home *** CD-ROM | disk | FTP | other *** search
/ Business Shareware / Business Shareware.iso / start / hobby / aa_51 / dms.c < prev    next >
C/C++ Source or Header  |  1993-02-09  |  10KB  |  551 lines
  1. /*  Utility programs for unit conversions and calendar
  2.  */
  3.  
  4. /*double PI = 3.14159265358979323846;*/
  5. #include "kep.h"
  6.  
  7. char *intfmt = "%d";
  8. char *lngfmt = "%ld";
  9. char *dblfmt = "%lf";
  10. char *strfmt = "%s";
  11.  
  12. /* Display Right Ascension and Declination
  13.  * from input equatorial rectangular unit vector.
  14.  * Output vector pol[] contains R.A., Dec., and radius.
  15.  */
  16. int showrd( msg, p, pol )
  17. char *msg;
  18. double p[], pol[];
  19. {
  20. double x, y, r;
  21. int i;
  22.  
  23. r = 0.0;
  24. for( i=0; i<3; i++ )
  25.     {
  26.     x = p[i];
  27.     r += x * x;
  28.     }
  29. r = sqrt(r);
  30.  
  31. x = zatan2( p[0], p[1] );
  32. pol[0] = x;
  33.  
  34. y = asin( p[2]/r );
  35. pol[1] = y;
  36.  
  37. pol[2] = r;
  38.  
  39. printf( "%s  R.A. ", msg );
  40. hms( x );
  41. printf( "Dec. " );
  42. dms( y );
  43. printf( "\n" );
  44. return(0);
  45. }
  46.  
  47.  
  48.  
  49.  
  50. /* Display magnitude of correction vector
  51.  * in arc seconds
  52.  */
  53. int showcor( strng, p, dp )
  54. char *strng;
  55. double p[], dp[];
  56. {
  57. double p1[3], dr, dd;
  58. int i;
  59.  
  60. if( prtflg == 0 )
  61.     return(0);
  62.  
  63. for( i=0; i<3; i++ )
  64.     p1[i] = p[i] + dp[i];
  65.  
  66. deltap( p, p1, &dr, &dd );
  67. printf( "%s dRA %.3lfs dDec %.2lf\"\n", strng, RTS*dr/15.0, RTS*dd );
  68. return(0);
  69. }
  70.  
  71.  
  72.  
  73. /* Radians to degrees, minutes, seconds
  74.  */
  75. int dms( x )
  76. double x;
  77. {
  78. double s;
  79. int d, m;
  80.  
  81. s = x * RTD;
  82. if( s < 0.0 )
  83.     {
  84.     printf( " -" );
  85.     s = -s;
  86.     }
  87. else
  88.     printf( "  " );
  89. d = s;
  90. s -= d;
  91. s *= 60;
  92. m = s;
  93. s -= m;
  94. s *= 60;
  95. printf( "%3dd %02d\' %05.2f\"  ", d, m, s );
  96. return(0);
  97. }
  98.  
  99.  
  100.  
  101.  
  102. /* Radians to hours, minutes, seconds
  103.  */
  104. #define RTOH (12.0/PI)
  105.  
  106. int hms( x )
  107. double x;
  108. {
  109. int h, m;
  110. double s;
  111.  
  112. s = x * RTOH;
  113. if( s < 0.0 )
  114.     s += 24.0;
  115. h = s;
  116. s -= h;
  117. s *= 60;
  118. m = s;
  119. s -= m;
  120. s *= 60;
  121. printf( "%3dh %02dm %06.3fs  ", h, m, s );
  122. return(0);
  123. }
  124.  
  125. /*        julian.c
  126.  *
  127.  * This program calculates Julian day number from calendar
  128.  * date, and the date and day of the week from Julian day.
  129.  * The Julian date is double precision floating point
  130.  * with the origin used by astronomers.  The calendar output
  131.  * converts fractions of a day into hours, minutes, and seconds.
  132.  * There is no year 0.  Enter B.C. years as negative; i.e.,
  133.  * 2 B.C. = -2.
  134.  *
  135.  * The approximate range of dates handled is 4713 B.C. to
  136.  * 54,078 A.D.  This should be adequate for most applications.
  137.  *
  138.  * B.C. dates are calculated by extending the Gregorian sequence
  139.  * of leap years and century years into the past.  This seems
  140.  * the only sensible definition, but I don't know if it is
  141.  * the official one.
  142.  *
  143.  * Note that the astronomical Julian day starts at noon on
  144.  * the previous calendar day.  Thus at midnight in the morning
  145.  * of the present calendar day the Julian date ends in .5;
  146.  * it rolls over to tomorrow at noon today.
  147.  *
  148.  * The month finding algorithm is attributed to Meeus.
  149.  *
  150.  * - Steve Moshier
  151.  */
  152.  
  153.  
  154. char *months[12] = {
  155. "January",
  156. "February",
  157. "March",
  158. "April",
  159. "May",
  160. "June",
  161. "July",
  162. "August",
  163. "September",
  164. "October",
  165. "November",
  166. "December"
  167. };
  168.  
  169. char *days[7] = {
  170. "Sunday",
  171. "Monday",
  172. "Tuesday",
  173. "Wednesday",
  174. "Thursday",
  175. "Friday",
  176. "Saturday"
  177. };
  178.  
  179. long cyear = 1986;
  180. extern long cyear;
  181. static int month = 1;
  182. static double day = 1.0;
  183. short yerend = 0;
  184. extern short yerend;
  185.  
  186. double getdate()
  187. {
  188. double J;
  189. double caltoj();
  190.  
  191. /* Get operator to type in a date.
  192.  */
  193. getnum( "Calendar date: Year", &cyear, lngfmt );
  194.  
  195. if( (cyear > 53994) || (cyear < -4713) )
  196.     {
  197.     printf( "Year out of range.\n" );
  198.     goto err;
  199.     }
  200. if( cyear == 0 )
  201.     {
  202.     printf( "There is no year 0.\n" );
  203. err:    J = 0.0;
  204.     goto pdate;
  205.     }
  206.  
  207. getnum( "Month (1-12)", &month, intfmt);
  208. getnum( "Day.fraction", &day, dblfmt );
  209.  
  210. /* Find the Julian day. */
  211. J = caltoj(cyear,month,day);
  212. /*printf( "Julian day %.1f\n", J );*/
  213.  
  214. pdate:
  215. /* Convert back to calendar date. */
  216. /* jtocal( J ); */
  217. return(J);
  218. }
  219.  
  220.  
  221.  
  222. /*     Calculate Julian day from Gregorian calendar date
  223.  */
  224. double caltoj( year, month, day )
  225. long year;
  226. int month;
  227. double day;
  228. {
  229. long y, a, b, c, e, m;
  230. int BC;
  231. double J;
  232.  
  233.  
  234. BC = 0;
  235. /* The origin should be chosen to be a century year
  236.  * that is also a leap year.  We pick 4801 B.C.
  237.  */
  238. y = year + 4800;
  239. if( year < 0 )
  240.     {
  241.     BC = 1;
  242.     y += 1;
  243.     }
  244.  
  245. /* The following magic arithmetic calculates a sequence
  246.  * whose successive terms differ by the correct number of
  247.  * days per calendar month.  It starts at 122 = March; January
  248.  * and February come after December.
  249.  */
  250. m = month;
  251. if( m <= 2 )
  252.     {
  253.     m += 12;
  254.     y -= 1;
  255.     }
  256. e = (306 * (m+1))/10;
  257.  
  258. a = y/100;    /* number of centuries */
  259. if( year <= 1582L )
  260.     {
  261.     if( year == 1582L )
  262.         {
  263.         if( month < 10 )
  264.             goto julius;
  265.         if( month > 10)
  266.             goto gregor;
  267.         if( day >= 15 )
  268.             goto gregor;
  269.         }
  270. julius:
  271.     printf( " Julian Calendar assumed.\n" );
  272.     b = -38;
  273.     }
  274. else
  275.     { /* -number of century years that are not leap years */
  276. gregor:
  277.     b = (a/4) - a;
  278.     }
  279.  
  280. c = (36525 * y)/100; /* Julian calendar years and leap years */
  281.  
  282. /* Add up these terms, plus offset from J 0 to 1 Jan 4801 B.C.
  283.  * Also fudge for the 122 days from the month algorithm.
  284.  */
  285. J = b + c + e + day - 32167.5;
  286. return( J );
  287. }
  288.  
  289.  
  290.  
  291.  
  292. /* Calculate month, day, and year from Julian date */
  293.  
  294. int jtocal( J )
  295. double J;
  296. {
  297. int month, day;
  298. long year, a, c, d, x, y, jd;
  299. int BC;
  300. double dd;
  301.  
  302. if( J < 1721425.5 ) /* January 1.0, 1 A.D. */
  303.     BC = 1;
  304. else
  305.     BC = 0;
  306.  
  307. jd = J + 0.5; /* round Julian date up to integer */
  308.  
  309. /* Find the number of Gregorian centuries
  310.  * since March 1, 4801 B.C.
  311.  */
  312. a = (100*jd + 3204500L)/3652425L;
  313.  
  314. /* Transform to Julian calendar by adding in Gregorian century years
  315.  * that are not leap years.
  316.  * Subtract 97 days to shift origin of JD to March 1.
  317.  * Add 122 days for magic arithmetic algorithm.
  318.  * Add four years to ensure the first leap year is detected.
  319.  */
  320. c = jd + 1486;
  321. if( jd >= 2299160.5 )
  322.     c += a - a/4;
  323. else
  324.     c += 38;
  325. /* Offset 122 days, which is where the magic arithmetic
  326.  * month formula sequence starts (March 1 = 4 * 30.6 = 122.4).
  327.  */
  328. d = (100*c - 12210L)/36525L;
  329. /* Days in that many whole Julian years */
  330. x = (36525L * d)/100L;
  331.  
  332. /* Find month and day. */
  333. y = ((c-x)*100L)/3061L;
  334. day = c - x - ((306L*y)/10L);
  335. month = y - 1;
  336. if( y > 13 )
  337.     month -= 12;
  338.  
  339. /* Get the year right. */
  340. year = d - 4715;
  341. if( month > 2 )
  342.     year -= 1;
  343.  
  344. /* Day of the week. */
  345. a = (jd + 1) % 7;
  346.  
  347. /* Fractional part of day. */
  348. dd = day + J - jd + 0.5;
  349.  
  350. /* post the year. */
  351. cyear = year;
  352.  
  353. if( BC )
  354.     {
  355.     year = -year + 1;
  356.     cyear = -year;
  357.     printf( "%ld B.C. ", year );
  358.     }
  359. else
  360.     printf( "%ld ", year );
  361.  
  362. day = dd;
  363. printf( "%s %d %s", months[month-1], day, days[a] );
  364.  
  365. /* Flag last or first day of year */
  366. if( ((month == 1) && (day == 1))
  367.     || ((month == 12) && (day == 31)) )
  368.     yerend = 1;
  369. else
  370.     yerend = 0;
  371.  
  372. /* Display fraction of calendar day
  373.  * as clock time.
  374.  */
  375. a = dd;
  376. dd = dd - a;
  377. hms( 2.0*PI*dd );
  378. if( J == TDT )
  379.     printf( "TDT\n" ); /* Indicate Terrestrial Dynamical Time */
  380. else if( J == UT )
  381.     printf( "UT\n" ); /* Universal Time */
  382. else
  383.     printf( "\n" );
  384. return(0);
  385. }
  386.  
  387.  
  388. /* Reduce x modulo 360 degrees
  389.  */
  390. double mod360(x)
  391. double x;
  392. {
  393. long k;
  394. double y;
  395.  
  396. k = x/360.0;
  397. y = x  -  k * 360.0;
  398. while( y < 0.0 )
  399.     y += 360.0;
  400. while( y > 360.0 )
  401.     y -= 360.0;
  402. return(y);
  403. }
  404.  
  405.  
  406.  
  407.  
  408. /* Reduce x modulo 2 pi
  409.  */
  410. #define TPI (2.0*PI)
  411. double modtp(x)
  412. double x;
  413. {
  414. double y;
  415.  
  416. y = floor( x/TPI );
  417. y = x - y * TPI;
  418. while( y < 0.0 )
  419.     y += TPI;
  420. while( y >= TPI )
  421.     y -= TPI;
  422. return(y);
  423. }
  424.  
  425.  
  426. /* Get operator to type in hours, minutes, and seconds
  427.  */
  428. static int hours = 0;
  429. static int minutes = 0;
  430. static double seconds = 0.0;
  431.  
  432. double gethms()
  433. {
  434. double t;
  435.  
  436. getnum( "Time: Hours", &hours, intfmt );
  437. getnum( "Minutes", &minutes, intfmt );
  438. getnum( "Seconds", &seconds, dblfmt );
  439. t = (3600.0*hours + 60.0*minutes + seconds)/86400.0;
  440. return(t);
  441. }
  442.  
  443. int getnum( msg, num, format )
  444. char *msg;
  445. double *num;
  446. char *format;
  447. {
  448. char s[40];
  449.  
  450. printf( "%s (", msg );
  451. if( format == strfmt )
  452.     printf( format, num );
  453. else if( format == dblfmt )
  454.     printf( format, *(double *)num );
  455. else if( format == intfmt )
  456.     printf( format, *(int *)num );
  457. else if( format == lngfmt )
  458.     printf( format, *(long *)num );
  459. else
  460.     printf( "Illegal input format\n"  );
  461. printf( ") ? ");
  462. gets(s);
  463. if( s[0] != '\0' )
  464.     sscanf( s, format, num );
  465. return(0);
  466. }
  467.  
  468.  
  469.  
  470. /*
  471.  * Convert change in rectangular coordinatates to change
  472.  * in right ascension and declination.
  473.  * For changes greater than about 0.1 degree, the
  474.  * coordinates are converted directly to R.A. and Dec.
  475.  * and the results subtracted.  For small changes,
  476.  * the change is calculated to first order by differentiating
  477.  *   tan(R.A.) = y/x
  478.  * to obtain
  479.  *    dR.A./cos**2(R.A.) = dy/x  -  y dx/x**2
  480.  * where
  481.  *    cos**2(R.A.)  =  1/(1 + (y/x)**2).
  482.  *
  483.  * The change in declination arcsin(z/R) is
  484.  *   d asin(u) = du/sqrt(1-u**2)
  485.  *   where u = z/R.
  486.  *
  487.  * p0 is the initial object - earth vector and
  488.  * p1 is the vector after motion or aberration.
  489.  *
  490.  */
  491.  
  492. int deltap( p0, p1, dr, dd )
  493. double p0[], p1[];
  494. double *dr, *dd;
  495. {
  496. double dp[3], A, B, P, Q, x, y, z;
  497. int i;
  498.  
  499. P = 0.0;
  500. Q = 0.0;
  501. z = 0.0;
  502. for( i=0; i<3; i++ )
  503.     {
  504.     x = p0[i];
  505.     y = p1[i];
  506.     P += x * x;
  507.     Q += y * y;
  508.     y = y - x;
  509.     dp[i] = y;
  510.     z += y*y;
  511.     }
  512.  
  513. A = sqrt(P);
  514. B = sqrt(Q);
  515.  
  516. if( (A < 1.e-7) || (B < 1.e-7) || (z/(P+Q)) > 5.e-7 )
  517.     {
  518.     P = zatan2( p0[0], p0[1] );
  519.     Q = zatan2( p1[0], p1[1] );
  520.     Q = Q - P;
  521.     while( Q < -PI )
  522.         Q += 2.0*PI;
  523.     while( Q > PI )
  524.         Q -= 2.0*PI;
  525.     *dr = Q;
  526.     P = asin( p0[2]/A );
  527.     Q = asin( p1[2]/B );
  528.     *dd = Q - P;
  529.     return(0);
  530.     }
  531.  
  532.  
  533. x = p0[0];
  534. y = p0[1];
  535. if( x == 0.0 )
  536.     {
  537.     *dr = 1.0e38;
  538.     }
  539. else
  540.     {
  541.     Q = y/x;
  542.     Q = (dp[1]  -  dp[0]*y/x)/(x * (1.0 + Q*Q));
  543.     *dr = Q;
  544.     }
  545.  
  546. x = p0[2]/A;
  547. P = sqrt( 1.0 - x*x );
  548. *dd = (p1[2]/B - x)/P;
  549. return(0);
  550. }
  551.